//Luokasta cNum peritty cRKN - Runge-Kutta-Nyström menetelmän toteutus.
//
//Joonas Koskinen
//Pro Gradu Työ
//Jyväskylän Yliopisto

#include "cRKN.h"

RKN::RKN() {
	double ivy, ivdy;
	setName("Runge-Kutta-Nyström");
	//tx and fx are cleared and set from cInfo when active potential is set.
	//If, for some reason, the active potential is not set and method is initiated we 
	//set tx and fx to 0 and 10 respectivly.
/*	ivy=0.0;
	ivdy=0.0;
	setStep(0.01);
//	set_tx(0.0);
//	set_fx(10.0);
	set_ty(ivy);
	set_fy(ivy);
	set_tdy(ivdy);
	set_fdy(ivdy);*/
}

RKN::RKN(double ivy, double ivdy) {
	setName("Runge-Kutta-Nyström");
	//tx and fx are cleared and set from cInfo when active potential is set.
	//If, for some reason, the active potential is not set and method is initiated we 
	//set tx and fx to 0 and 10 respectivly.
	setStep(0.01);
//	set_tx(0);
//	set_fx(10);
	set_ty(ivy);
	set_fy(ivy);
	set_tdy(ivdy);
	set_fdy(ivdy);
}

void RKN::methodStep(double& step, vector<double> &xx, vector<double> &yy, vector<double> &dy, cDY *DY, cPot *V) {
	double k1,k2,k3,k4;
	k1 = 0.5*step*DY->DE(xx.back(),V->potential(xx.back()),yy.back(),dy.back());
	k2 = 0.5*step*DY->DE(xx.back()+0.5*step,V->potential(xx.back()), yy.back()+0.5*step*(dy.back()+0.5*k1), dy.back()+k1);
	k3 = 0.5*step*DY->DE(xx.back()+0.5*step,V->potential(xx.back()), yy.back()+0.5*step*(dy.back()+0.5*k2), dy.back()+k2);
	k4 = 0.5*step*DY->DE(xx.back()+step,V->potential(xx.back()), yy.back()+step*(dy.back()+k3),yy.back()+2.0*k3);

	yy.push_back(yy.back() + step * (dy.back() + ((k1 + k2 + k3) / 3.0)));
	dy.push_back(dy.back() + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 3.0);
	xx.push_back(xx.back() + step);

}

RKN::~RKN(){}
